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SUMMARY 

The vi scous/invi scid interaction over transonic airfoils with and without 
suction is studied. The streamline angle at the edge of the boundary layer is 
used to couple the viscous and inviscid flows. The potential flow equations 
are solved for the inviscid flow field. In the shock region, the Euler 
equations are solved using the method of integral relations. For this, the 
potential flow solution is used as the initial and boundary conditions. An 
integral method is used to solve the laminar boundary-layer equations. Since 
both methods are integral methods, a continuous interaction is allowed between 
the outer inviscid flow region and the inner viscous flow region. 

To avoid the Goldstein singularity near the separation point the laminar 
boundary-layer equations are derived in an inverse form to obtain solutions 
for the flows with small separations. The displacement thickness distribution 
is specified instead of the usual pressure distribution to solve the boundary- 
layer equations. The Euler equations are solved for the inviscid flow using 
the finite volume technique and the coupling is achieved by a surface 
transpiration model. A method is developed to apply a minimum amount of 
suction that is required to have an attached flow on the airfoil. The suction 
parameter is varied based on the velocity profile parameter and the suction 
distribution obtained is considered to be close to the optimum value. The 
turbulent boundary layer equations are derived using the bi-logarithmic wall 
law for mass transfer. The solution method is similar to the laminar inverse 
boundary-layer approach. The results are found to be in good agreement with 
available experimental data and with the results of other computational 
methods. 
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Chapter 1 


INTRODUCTION 

Transonic flows are the flows in which the local flow speed is 
close to the local speed of sound. These flows occur m nozzles, over 
propellers and turbine blades, around blunt bodies flying supersonically 
and near airplanes which fly close to the Mach number of one. The 
interest in transonic flow started due to the problems encountered in 
the attempts to design efficient commercial aircraft which fly close to 
but below speed of sound. 

The most distinguishing feature of transonic flows is their mixed 
flow character. The acceleration of the initially subsonic flow over 
the forward portion of an airfoil is sufficient to provide an embedded 
region of supersonic flow adjacent to the airfoil surface. This 
supersonic region is terminated by a shock wave that recompresses the 
flow. 

The qualitative behavior of lift and drag coefficients (C^ and 

* 

C.) as functions of free stream Mach number M is discussed m [1J. 
a cd 

The critical Mach number M is the value of M for which an 

cr * 

embedded supersonic region first appears. As the Mach number increases 


* 


The numbers in brackets indicate references. 
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beyond M cr » the supersonic region grows, increasing the strength and 

extent ot the terminating shock; C also increases and c, essentially 

t a 

i t 

remains constant. As increases beyond M^, the drag rise Mach 

number, shock and viscous influences cause a rapid increase in drag and, 

eventually, a decrease in lift. Therefore, the optimum cruise Mach 

number is for the value of M just above M.. 

m J d 

The main objectives for fighter type aircraft are high lift at low 
drag level, high thrust-to-drag ratio for acceleration, and high load 
factors for maneuvers. These features make the analysis of transonic 
flow fields one of the most studied in fluid dynamics. To improve these 
factors limiting the performance of aircraft, detailed studies 

(comprising both wind tunnel testing and fluid dynamic computations) 
have to be performed. The high cost of transonic wind tunnel test time 
severely limits the number of configurations that can be considered in 
the search for the optimum design. Considerable attention has been 
directed m recent years in developing theoretical models for the design 
of aerodynamic bodies to supplement wind tunnel simulations. Theoreti- 
cal methods reduce the cost of wind tunnel testing and also provide a 
data base to check errors due to wind tunnel wall interference, in- 
correct Reynolds number scaling, model and sting deflections under 
load. A natural approach to the computation of transonic viscous- 
inviscid interaction is to use the Reynolds averaged Navier-Stokes 
equations as a global solution procedure. However according to present 
trends ( 2 ) , it could be a long time before computers will attain the 
power required to perform the routine engineering calculations using the 
Navier-Stokes equations. An alternative is to use a zonal solution 
method to compute the transonic viscous-mviscid interaction. This 


i 
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method involves solving different numerical algorithms tor various 
regions of the flow. An advantage of zonal methods is that simpler 
equations are used in regions where permissible and, consequently, 
computer requirements are reduced. A disadvantage of zonal methods is 
that information must be exchanged at zonal boundaries or zonal overlap 
regions, which can cause convergence and stability problems. 

The first major breakthrough to solve the governing nonlinear 
partial differential equations was made by Magnus and Yoshihara [3] . 
They presented a method to solve the steady supercritical planar flows 
over lifting airfoils using the finite-difference technique. Steady- 
state solutions were obtained as the asymptotic flow for large times. 
Mathematically, the description of steady transonic flows requires the 
solution of "mixed equations," which are elliptic in subsonic regions 
and hyperbolic in supersonic regions. The problem is nonlinear, and 
solutions generally contain discontinuities representing shock waves. 
Murman and Cole 1 1 3 presented a governing steady transonic small 
disturbance potential equation using a mixed finite-difference system. 
They used central difference formulas in the subsonic region, where the 
governing equation is elliptic and upwind difference formulas in the 
supersonic region, where it is hyperbolic. Krupp and Murman [4] 
extended this method to lifting airfoils. These initial successes 
started a new field of study, the computational transonic aerodynamics. 

After the approach suggested by Murman and Cole 11] for solution to 
the transonic flow problems there was considerable research interest to 
solve potential flow equations using different techniques. Lomax et al. 
[5], Klunker and Neuman 16], Schmidt et al. 17], Albone et al. [8], and 
Van der Vooren et al. 19] have used various forms of the transonic small 



4 


disturbance equations with the aim of improving the solutions. There is 
also a parallel effort to solve the exact potential equation. 
.Garabedian et al. [10] solved the full equation for velocity potential 
using a relaxation technique similar to Murman's transonic small 
perturbation method. In this method the central difference scheme for 
subsonic flow is replaced by a 'retarded' differencing scheme for 
supersonic flow. Jameson [11], Arlinger [12], Baker [13], and Caughey 
and Jameson [14] have obtained results for two-dimensional and axially 
symmetric flows. Jameson [15] solved the transonic full potential flow 
equation using a 'rotated* differencing scheme to confirm with the local 
stream direction. Carlson [16] solved the full, inviscid perturbation 
potential equation m a Cartesian system. The second-order partial 
differential equation is replaced with a nonconservative system of 
finite difference equations which includes Jameson's rotated 
differencing scheme at the supersonic points. Using a stream function 
approach, Hafez and Lovell [17] solved the transonic inviscid flow 
equations. Also numerous techniques were developed to improve the 
efficiency of the basic algorithms. Ballhaus et al. [18] developed 
implicit approximate factorization algorithms for steady-state transonic 
flow problems to accelerate convergence. South and Brandt [19] used a 
multigrid method to achieve this goal. A hybrid Poisson solver/ 
successive line over relaxation (SLOR) scheme proposed by Jameson [20] 
substantially increased the computational efficiency over the 
conventional SLOR scheme. 

The comparison between the solutions using the potential flow 
equations and experimental data is not very good in the mid to upper 
transonic regions. This is because of the isentropic assumption 
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inherent in these formulations. Recently, it was also discovered that 
these solutions are not unique; for certain angles of attack and Mach 
numbers, the full potential flow equation yields multiple solutions. 
Also, large negative or positive lift is predicted even for symmetric 
airfoils at zero angle of attack. Steinhoff and Jameson [21] noted that 
the nonuniqueness is not because of the discretization error, but 
because of the governing partial differential equation itself. Several 
investigators studied the occurrence and behavior of this anomaly, and 
concluded that it shows up in the conservative formulation of the 
partial differential equation 122-24], but not m the nonconservative 
formulation. This suggests that the anomaly is associated with the 
approximate treatment of shock waves within the potential formulation 
since this is the main difference between the conservative and 
nonconservative formulations. 

The solution of Euler equations do not show the anomalous behavior 
as exhibited by the potential solutions. However, the cost of computa- 
tions is higher. For flows in the upper transonic range it is necessary 
to solve the Euler equations because of the rotational nature of the 
flow. The comparison between the Euler and potential solutions is found 
to be excellent in the lower transonic or subsonic range. Using Euler 
equations in the whole flow field or in the shock region should improve 
the results. Euler equations can be solved either by using the method 
of integral relations (MIR) or by the finite-difference technique. The 
method of integral relations has been used by Holt and Masson [25], 
Melnik and Ives [26], Sato [27), and Tai [28] to solve the Euler 
equations. This method is valid for transonic flows with moderate to 

strong shocks. The disadvantage of MIR is that the solution procedure 
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requires man-machine interactions because of the multiple iterative 
processes involved. Liepmann {29] and Ackeret et al. (30) investigated 
the phenomenon of the transonic viscous-mviscid interaction. Their 
experiments showed that the shock and the boundary layer interact 
strongly with each other. This phenomenon is of great complexity 
because the behavior of the boundary layer depends mainly on the 
Reynolds number, whereas, the conditions in a wave are primarily 
dependent on the Mach number. The pressure disturbances caused by the 
shock propagate upstream through the subsonic portion of the boundary 
layer causing the flow to separate ahead of the shock. 

Bauer et al. 131 ,32] incorporated the Nash and Macdonald [33] 
turbulent boundary-layer method into the inviscid Garabedian method 
[10]. The viscous-inviscid interaction was taken into account using a 
solid displacement model. Collyer and Lock [34-36] used the lag 

entrainment method of Green et al. [37] to calculate the turbulent 

boundary layer. The surface transpiration model was used to represent 

the displacement effect of the boundary layer and wakes on the 
equivalent inviscid flow. This method has an advantage over the 
previous method in that the computational grid needs to be generated 

only once. Melnik [38] used a 'multi deck' model near the trailmg-edge 
region, based on the asymptotic theory of turbulent shear flows in the 
limit of infinite Reynolds number. The matching between the inviscid 
and viscous solutions is achieved using the surface transpiration model. 

Klineberg and Steger [39] treated the viscous-inviscid interaction 
by using a boundary- layer integral approach combined with a finite- 
difference relaxation technique for the small disturbance equations. 
Inviscid and viscous flows are treated separately even for strong 
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interactions. Brilliant and Adamson (40) used the method of matched 

asymptotic expansions for an incident shock interacting with an 
unseparated laminar boundary layer in transonic flow. Tai 141] coupled 
the inviscid transonic solution obtained by using the method of integral 
relations with the integral method developed originally by Lees and 
Reeves [42] and refined by Klineberg and Lees [43] for compressible, 

attached and separated laminar boundary layers. 

All of the methods discussed above are direct methods where the 

external pressure distribution is prescribed and the boundary-layer 
quantities are calculated. These methods exhibit the Goldstein 
singularity near the separation and are limited to attached flow 

conditions. As demonstrated by Goldstein [44], the boundary layer 
growth rate becomes infinite when the shear stress gradient is infinite 
causing a singularity at the point of separation. This was confirmed by 
Klineberg and Steger [45], Werle and Davis [46], and Pletcher and Dancey 
[47]. It is also posssible to get nonunique solutions from direct 
boundary- layer solutions. For non-similar flows this problem was noted 
by Murphy and King [48]. Catherall and Mangier [49] pointed out that 
this does not limit the validity of Prandtl's boundary-layer equations 
past separation. Using an inverse approach where the boundary- layer 
thickness distribution is specified, this problem could be avoided. 
Many of the recent researchers have used inverse boundary- layer 
equations to calculate separated flows. Carter [50] showed that the 
inverse boundary- layer solutions compare well with the Navier-Stokes 
solution. Cebeci et al. [51] used a nonlinear eigenfunction formula- 
tion, and Klineberg et al. [45] and Horton [52] specified shear stress 
in their calculations to obtain regular solutions. 
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The high cost and time requirements and the resources available for 
conducting experiments either m the wind tunnel or in flight warrant 
. the need for an approximate suction distribution from a reliable 

numerical method. Also, the solution procedure should be faster and 
efficient to perform the routine calculations on any desired airfoil. 
It is essential to have an idea about the suction quantities or the 

velocities to maintain full chord laminar flow to reduce viscous drag. 
This will lead to a better design of Natural Laminar Flow (NLF) or 

Laminar Flow Controlled (LFC) airfoils. The computational procedure 
should include both the laminar and turbulent boundary- layer models with 
a transition criteria. In the viscous-inviscid interaction near the 
shock and the trailing edge, the wake curvature effect have to be 
considered. The specific objectives of the present study, therefore, 
are to develop a computational method which would consider laminar as 
well as turbulent attached boundary-layer interaction for the flow over 
transonic airfoils, and to develop a method to obtain the suction 

distribution for maintaining attached flow on the airfoil. 

In Chap. 2 governing equations for the inviscid flow as well as the 
direct and inverse boundary layer equations are presented. Method of 
solution to solve these equations is discussed in Chap. 3. Also, the 
interaction models to couple the outer inviscid flow with the inner 
viscous flow are presented. Results are obtained for several airfoils 
using the direct and the inverse boundary- layer approaches for laminar 
and turbulent flows. These are compared with the available experimental 


as well as other numerical results in Chap. 4. 



Chapter 2 


THEORETICAL FORMULATION 


2.1 Basic Formulation for Inviscid/Viscous Interaction 

The mixed flow field over an airfoil in the transonic range is 
illustrated in Fig. 2.1. Except in the shock region the potential flow 
equation is solved by a finite-difference relaxation technique. Euler 
equations are solved in the shock region. In Sec. 2.1.1 the governing 
equations for the mviscid flow are presented. The direct boundary- 
layer equations are derived in Sec. 2.1 .2 and the corresponding initial 
conditions are given in Sec. 2.1.3. 

2.1.1 Inviscid Flow 

Under the assumption of mviscid and irrotational flow the 
transonic potential flow equation is given by [53] 

fa 2 - - 26 6 d> + (a 2 - 6 2- )* =0 (2.1) 

v x ' xx x y xy v y ' yy 

where $ is the velocity potential. 

Introducing a perturbation potential of the form 
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Equation (2.1) can be written as 


(a 


2 

u )« 


XX 


- 2uv$ + (a 
xy 


v 2 )$ 


yy 


= o 


(2.3) 


where u 


9 and v 
x 


♦ are the velocity components, and 


2 

a 




2 

v 



(2.4) 


Equation (2.3) has the form locally, either of a wave equation (hyper- 
bolic type) representing supersonic flow ( 9 ^ > 1 ) , or of a Laplace 
equation (elliptic type) representing the subsonic flow ( <J>^ < 1). The 

nonlinear term, uv$ allows the transition from one type to another. 

xy 

The boundary conditions for mviscid flow are given as 



v u •'body 


(2.5a) 


r — i 

♦ = - — tan f B tan( 0 - a) ) 

2 7T 


(2.5b) 


where circulation T is determined by the change in potential across 
the Kutta-Joukowski cut at the trailing edge (TE), i.e.. 


T * (♦ -9 ) (2.5c) 

y=0 y=0~ TE 

The potential flow equation is rearranged in rotated coordinates 
parallel and perpendicular to the local velocity. This rearrangement 
permits coordinate stretching m the physical plane and avoids computa- 
tional problems in the supersonic region. Several methods are available 
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[32,54] to solve the potential Eq. (2.3) subjected to the boundary 
condition, Eq. (2.5). 

, The isentropic assumption inherent in the potential flow equation, 
Eq. (2.3), is not valid m the case of moderate to strong shocks across 
which the increase in entropy cannot be neglected. Therefore, in the 
shock region it is necessary to solve the Euler equations for the 
inviscid flow solution. The Euler equations are expressed in vector 
form as 


where 


A + B = 0 
x y 



- pu - 


- pv 

— 

2 



A » 

pu + p 

; B = 

puv 


_ puv . 


Lpv + pj 


and 


2 2 

C T + Vo (u + v ) = constant 
P z 



(2.6a) 


(2.6b) 


(2.6c) 


2.1.2 Viscous Flow 

The governing equations for a steady, two-dimensional compressible 
laminar boundary layer in coordinates parallel and normal to the surface 


are 



1 3 


Continuity: 


3( pu) 3( pv) 

— + — = o 

3s 3n 


(2.7) 


s-momentum: 


3u 3v dp 1 3 r 3u > 

Ts + *" = “ K is + iT UT (•* to) 


( 2 . 8 ) 


n-momentum : 


3p 

-f- = 0 
3n 


(2.9) 


Energy: 


2 2 

C T + V? (u + v ) = 0 
P A 


( 2 . 10 ) 


Equation of state: 


P = PRT 


( 2.11 ) 


Illingworth [55] and Stewartson [56] have shown that for an 
adiabatic surface and for a fluid with a Prandtl number of unity, the 
compressible boundary- layer equation can be transformed into an 
incompressible form. By applying the Stewartson transformation, 


s a P 
e e 


6- / f f* s> n- / 

0 00 CO 


n a p 


e _e j>_ 

0 a oo p oo p e 


dn 


(2.11a) 


P a _ , n p a 

U = ( a /a )uj V = f— — — ) u "T f — — — •&— dn 
»' e V P M a ' 3s ' p a p 

e e 0 “ » e 


P a 

Co €0 P 

Pap 
e e e 


(2.11b) 
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Eqs. (2.7) - (2.9) are reduced to the incompressible form as 


3U j)V 

H + 


0 


( 2 . 12 ) 




dU 


1 3 U 


Re 


3n 


(2.13) 


In order to avoid the semi-empirical features inherent in the methods 
such as Crocco-Lees (57) , a moment of momentum equation is used in 
addition to Eqs. (2.12) and (2.13). Upon multiplying Eq. (2.13) by U, 
one finds 


2 3U 3U dU e 

U H + W IS = UU e df 


Re 


_3 U 

3n 2 


(2.14) 


The governing partial differential equations, Eqs. (2.12) - (2.14), 
are integrated across the boundary layer resulting in three ordinary 
differential equations. These equations can be written in the matrix 
form as 


* 

6 . 

i 


* dJ 
{ i dH 


* 

6 . 

1 A 
M 


6. (2H + 1) 
M 


3 6* J 

l 

M 


* 

*d6.“ 
i 

ds 


dH 

ds 


dM 
€ 

ds 


JL 

m 


1 + m 


1 + m 


M ^ V 

b — 2 — + 

P M U 

e Re t e 

6 

l 


M 

00 

6 r 


u 


e Re 


w 


tan G - 


(2.15) 
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where 


e 

Hr. c 
« 

i 


— f— ) 

u ' 3n * 


3n -'rpO ' 


H = 


26 


i. t r 3U \ 2 , 

ZT [ dT1 


1 + m 


F = H + 


m 


2=4/ 

6 O U e 
l 


Re 


V „ * 

= — c 6 
v i 


(2 ♦ j«4 

1 Vi 


m 


1+m 


3~Y~1 

T-1 


-)h + •= : -^-r + 


2 

M - 1 
e 

m (1+m ) 
e e 


For a given velocity profile — , the integrals in R and Z can be 

e 

evaluated and the system of Eq. (2.15) can be solved simultaneously. 

Lees and Reeves 142) have shown that the solutions of the Falkner- 
Skan [58) equation for similar flows including reversed flow profiles 
calculated by Stewartson [591 can be used to determine J, Q, H, R, 
and Z as functions of a single parameter ‘a’. This is referred to as 
the velocity profile parameter and is given by 


a 


3(U/U e ) 


0.99 


f"(0) 


0 < a < 4 

for attached profile 



I b 


a 



e 


“>£•. 0 


0.99 


0 < a < 1 

for separated profile 


a 


(-) 

e dividing 
streamline 


= (f * ) 


f =0 ' 


0 < a < 0.46 

for wake reversed profile 


a 



n=o 


(f 1 ) 


c=0 ' 


0 < a < 1 

for wake forward profile 


The quantity f is obtained from the solution of the Falkner-Skan 
equation 

f ... + ff « + b(1 _ f . 2 j = o (2.16) 


subjected to the boundary conditions 

f(0) = 0 

f'(0) = 0 

lim f * ( = 1 


These relations 

for 'a ' 

are valid 

not only 

in the 

case 

where 

there is 

no mass 

transfer at the 

wall. 

but 

also m 

cases 

with 

mass 

transfer 

provided 

f 

is obtained by 

solving Eq. 

(2.16) 

with 

the 

boundary 


conditions 

f(0) - E 

f'(0) = 0 
lim 


f'(-) = 1 
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Considering e to be a small perturbation, Eq. (2.16) can be written as 

(f 1 ♦ ef 2 > ' ' ' + (f 1 + + ef 2 >" + b(l - (f 1 + ef 2 ) * 2 ) = 0 (2.17) 

and the boundary conditions as 

f^O) + ef 2 (0) = e (2.17a) 

f'(0) + ef 2 (0) - 0 (2.17b) 

lim (f'U) + Ef^(0)= 1 (2.17c) 

0 

For small perturbations, the e equation is expressed as 

f"* + f^’ + b(l - f’ 2 ) * 0 (2.18) 

f^O) - 0 
f^(0) « 0 

and 

lim f 1 ( t) «= 1 
C+* 1 

The t equation can be expressed 

f”’ + f 2 f 1 + f 1 f 2 + 2bf 1 f 2 * 0 (2.19) 


and the boundary conditions are 



IB 


f (0) = 1 
2 

f^(0) = 0 
lim f *(;) = 0 

C-h» 2 

Equation (2.18) is the Falkner-Skan equation and Eq. (2.19) is the 
auxiliary equation governing {^. 

Solutions to the Falkner-Skan equation are available in the 
literature for different pressure gradients. But to evaluate the 
suction quantities, very accurate solutions to the Falkner-Skan equation 
and the auxiliary equation due to suction are needed. Therefore, these 
equations subjected to the respective boundary conditions are solved by 
using the state variable approach of Forbnch [60]. A solution accurate 
to sixth decimal place can be obtained using this method. 

The integral quantities used in Eq. (2.15) consists of the 
perturbed and unperturbed parts. The unperturbed integral quantities 
are, 

C 0.99 

H " * / f .'Cl ~ (2.20a) 

6 0 
il 

i ?0 - 99 , 

j, - / f ;0 - f; )n (2.20b) 

6 ii 0 

• 

Q ■ 6. f" (0) 

W 1 il \ 


(2.20c) 
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* 0*99 

R , - 26 u / f ; «•« 


jr IVW - f i (0) l 


11 


where 


* 

6 i1 


^0 .99 

/ 0 - *;) <n 


The perturbed integral quantities are 


«2 =~ 


6 i (1 + b) 


[(1 + b) ( 1 + H 1 )(f 2 (O 0>99 - 1) 
+ f 2 (0) - f 2 ( ^0.99l 


0.99 


J 2 -V I" ♦ J ,>(£ 2 U) 0-99 - 0 - 3 / £J 2 £I,d C } 


e 2 - - £;■ (o) (f 2 (c ) 0-99 - i) 


0 . 99 


R 2 --i- (' - f 2 U> 0 . 99 ) + 4«* J £-f" 2 dt 


fi. 

1 


(1 + z. ) 


* [^2 ^ ^ a . Q<3 “ 1 ] 


'2 * 1 2"* 0.99 

0. 

1 


All these quantities are expressed as polynomials in 'a.' 
combination of the velocity profile parameter and the suction 


( 2 . 20d ) 
(2.20e) 


( 2 . 20f ) 

(2.20g) 

(2.20h) 

(2.20i) 

(2.20j> 

For any 
parameter. 
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the total quantities can be evaluated as H = + eH^ , J 

etc. 


J 1 + eJ 2‘ 


Nov the system of Eqs. (2.15) can be expressed as 


d6* 

1 

= D 1 

d 

1 

3 1 2 

a i3 

l 

ds 

d 2 

a 22 

a 23 



d 3 

a 32 

a 33 


( 2.21 ) 




3 U 

d 1 

*13 

da 

1 


J 


ds 

_ dH 

a 21 

d 2 

*23 


1 da 

a 31 

d 3 

a 33 


dM 
e 

ds 


a 

a „ 

1 1 

12 

a 

a _ 

21 

22 

a 

a 

31 

32 



( 2 . 22 ) 


(2.23) 


where a... d and D are given m Appendix A. 

ID 3 1 

In the usual boundary-layer calculations, the pressure distribution 
along the airfoil is computed using an inviscid method, and the 
boundary- layer quantities are evaluated. Since the solution of the 

viscous and inviscid equations is not simultaneous, mass transfer 
between these two regions is not allowed and only the momentum and 
moment of momentum equations are used. 
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For a specified pressure distribution, Eq. (2.15) reduces to 


dS 

i 

ds 


M „ V 6 (2H + 1) dM 

s _f! £L + _i e 

P M * U m ds 

e Re g e e 

l 


(2.24) 


** dJ, 

6 1 dH 1 


dH 

ds 


M 

o_ 

6 ~ 


3 6 J dM 
l e 


M Re * 
6 6 

l 


M 


ds 


This may be rewritten as 




where 


da 

ds 


1 


D 

2 da 



(2.25) 


(2.26) 


D 


2 


a 


a 


21 

S 22 

31 

a 32 


b 


2 


dM 

d 2 " a 23 ds~ 


b 


3 



dM 
e 

a 33 ds 


These equations have a singularity at the separation point. If suction 
is not applied before the flow separates, the denominator becomes zero 
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and the boundary-layer equations canndt be applied past separation 
point. This formulation, termed as the weak interaction formulation, 
.should be applied only in the attached flow region. 


2.1.3 Initial Conditions 

To avoid the stagnation point singularity the calculations are 
started away from the leading edge. Klmeberg and Steger [61 J, Tai 
[28], and Ram et al. 162] assumed that the flow is locally similar to 
derive the initial conditions. Although this assumption is valid for 
only thin airfoils, it was found by Tai [28] that these conditions can 
be applied also to blunt airfoils to obtain converged solutions. 

For a locally similar flow, 4^- = 0 and Eq. (2.26) reduces to 

ds 



b 


2 


0 


b 


3 


or 


a b - b a 
21 3 2 31 


0 


or 


3 6 J dH 
i e 

M ds 
e 


(2.27) 


dM 

Upon substituting for . ^ from Eq. (2.27) into Eq. (2.25), there is 

ds 

obtained 


M V M V 

® R 0 ) •» 0 0 ) 

H % „ . + U =J8 M + U 

e Re * e e Re * e 

6 i 6 i 


6. ( 2H+1 ) dM 
1 e 


M 


ds 
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R + 2 HR - 3J0 
J(1 - H)Re 5 * 


(1 + H - 3 J ) 
J(1 - H) 


V 

0) 


6 

e 


(2.27a) 


or 


* 

6 


l 



K, 6. + K_ 
1 l 2 


(2.27b) 


where 


= (1 + H - 3 J ) V o) 
K 1 " J(1 - H) U 

e 


_ (R 4- 2 HR - 3JQ) 
K 2 J(1 - H)Re 


An integration of Eq. (2.27b) with respect to S gives 


* 



or 


[S* K, - Kj «n(K, . K 2 )] 0 - IK, S1 0 


(2.28) 


From Eq. (2.27), one finds 


JQ - HR 1 

J(1 - H)Re dM 
c e 


v 2 


ds 

Starting with an initial value of 'a', the velocity profile quantities 

are computed and both sides of Eq. (2.28) are evaluated. Then the 

initial value of 'a' is updated depending upon the ratio of the right- 

hand side to the left-hand side. This iteration process is continued 

until there is no appreciable change in the velocity profile parameter 

* 

'a'. The initial displacement thickness 6^ is calculated for this 


value of 
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2*2 Laminar Inverse Boundary-Layer Interactions 

, The transformed Euler equations are presented in Sec. 2.2.1 to 
solve for the entire inviscid flow. In Sec. 2.2.2 the boundary- layer 
equations are written in an inverse form to avoid the Goldstein 
singularity past the separation point and the corresponding initial 
conditions are derived in Sec. 2.2.3. 

2.2.1 Inviscid Flow 

In the direct boundary-layer computations, Euler equations are 
solved only in the shock and wake regions. Because of the multiple 
interative process involved in solving the two different types of 
equations for the inviscid flow solution, the computational procedure is 
difficult. It is convenient to solve the Euler equations in the entire 
inviscid flow region especially when the weak interaction equations are 
solved in the boundary layer. 

For numerical calculations Eqs. (2.6) are transformed into general 
coordinates using the transformation 

R = R(x,y) 

S = S (x,y) 

This results in 

A + B «= 0 (2.29) 

R S 


where 



A = (R A + R B)/J 
R x y c 


B = (S A + S B)/J 
S x y c 


J = (R S - R S ) 
c x y y x 


2.2.2 Inverse Boundary-Layer Equations 

In the usual boundary-layer method the boundary- layer quantities 
are calculated for a specified pressure distribution. These methods, 
termed as direct methods, exhibit the Goldstein singularity near the 
separation. The flow past separation cannot be calculated because of 
this singularity. However, this behavior does not limit the validity of 
the Prandtl equations. In an inverse boundary-layer method the 
boundary-layer thickness is specified and the pressure distribution is 
evaluated. An integral method is used to solve the inverse boundary- 
layer equations in the present approach. For a known boundary- layer 
thickness distribution, the viscous governing equations, Eqs. (2.15), 
reduce to 



* dJ 
*1 dH 


• 


— — 


«• — . 

* 




* 

6 (2H + 1) 




M _ V d 6 

1 


dH 


a " Q . » H h 

M 


ds 


0 M Re + U H ds 

e 




e e 





i 

* 




* 

36. J 


dM 


M _ V d 6 

i 


e 


3 

8 

% 

le 

4 

|h 

M 


ds 


M Re . U ° ds 

e 




e 6 6 

_ 


• m 


1 _ 


This is written in a compact form as 



26 


a 22 

*23 


dH " 
ds 

S 

K 3 




dM 



*32 

*33_ 


6 

ds 


K 4_ 


d«* d 6* 

where K. = d - H — — and K = d_ - J — " . From Eq. (2 
3 2 ds 4 3 ds M 

follows that 


da 

ds 


1 

_ dH 

D • 

3 da 



dM 
e 

ds 


1 



22 


‘32 



where D *= a a - a a . Equations (2.31) and (2.32) 

■3 1 4 J J Jb 

expressed in an alternate form as 


da 

ds 


dH 

da 


BM 6. 


M 2 Re 


(3JQ - 


* * 

6 d 6 

R(2H + 1)) + (J<2H + 1) 

e 


where 


dM 
e 

ds 


JL 

D. 


BM 6. 
e l 

M Re ^ 


dS 




*2 


(3J - (2H + 1) 


M 


dH 


(2.30b) 

30b) it 
(2.31 ) 

(2.32) 
can be 

• 3JH ) 

(2.33) 

(2.34) 
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The 

inverse boundary-layer 

equations, Eqs. 

(2. 

33) 

and (2.34) 

, are 

checked 

for 

singularities for 

both the attached 

and 

separated 

flow 

conditions. 

However, for extensive regions 

of 

separation, 

these 


equations are not valid. 


2.2.3 Initial Conditions 

Assuming the flow to be locally similar near the leading edge, 
da 

i.e., — — = 0, Eq. (2.33) reduces to 
ds 


or 


where 


M Re 

6 6 * 
l 


* 

d 6 

(3JQ - R( 2H + 1 ) ) + (j(1 


H) ) = 0 


d 6* ^ 

ds = M Re K 5 /K 6 
e ,* 


(2.35) 


K 5 = 3JQ - R(2H + 1) 


K = J ( H - 1 ) 

D 

* 

d 6 

Upon substituting for - in Eq. (2.34), there is obtained 



where 


(2.36) 


D . - 3J 

4 


(2H + 


1 ) &- 
’ dH 
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Since the velocity profile parameter and the profile quantities are 
locally constant, Eq. (2.36) can be integrated to obtain 


M 

e 



/ „ dJ> 1 

^ ® dH ^ + * 


&M 


* Re 
6 * 

l 6. 

i 


v dH 



(2.37) 


A substitution for M e in Eq. (2.35), results in 

* 

d6. 

3T 7 - c mK + (" m - K s d 4 ‘ 2 - 38 > 

0 , 

1 

Solving for the velocity profile parameter 'a', which satisfies Eq. 
(2.38), gives the initial value for 'a' and for that value the initial 
Mach number is calculated from Eq. (2.35). 


2.3 Turbulent Inverse Boundary-Layer Interaction 


The governing equations for a compressible, turbulent boundary 
layer in coordinates parallel and normal to the surface are 163] 


Continuity: 

3( pu) 3( pv) 
3s 3n 


(2.39) 


Momentum: 




, dp 1 3 r 3u \ 

k ds + p V c 3n 3n^ 

00 00 


(2.40) 
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By applying the Stewartson transformation Eqs. (2.39) and (2.40) are 
reduced to the incompressible form as 


Continuity: 


U + V = 0 

C n 


(2.41 ) 


Momentum : 


UU + VU 
C h 


U (U ) + — 

e e z Re 


(B 1 U n } n 

i n n 


(2.42) 


The input quantities are transformed to the incompressible form and 
after the boundary-layer calculations, the results are transformed back 
into the compressible plane. 


2.3.1 Inverse Boundary-Layer Integral Method 

The governing partial differential equations are integrated across 
the boundary layer as 


6 

/ u dn + 

0 c 


6 

/ v dn 
A n 


o 


(2.43a) 


or 



dn 


v 

u 


V 

e 


6 

/ K 
0 


VU - U (U ) - 

n e e c 


1 3 r a 8U \ -i , 

ST ( 8 i 


(2.43b) 


(2.44) 


A substitution for V from Eq. (2.41) in Eq. (2.44) gives 






(2.45) 



3(J 


As in the laminar boundary-layer problem, a moment of momentum equation 
is used to obtain the closure relationship. 


Velocity Profile 

The velocity profile expression similar to Kuhn and Nielson [64] 
has been used to eliminate the n dependence of the integral equations 
and this is given by 


= 2.5 £n( 1 + n + ) + 5.1 - (3.39n + + 5.1)e"°' 37n 

T 

(2.46) 

+ \ " cos(n *)J + vj2.5 ln( 1 + n + ) + 5.1 f 


The parameter 


U 

T 


is the friction velocity and is given by 


U = 
T 


•VlTj’llTj/rt* 


(2.47) 


Equation (2.46) consists of an inner part, consisting of a laminar 
sublayer and the law of the wall function, and an outer part, a wake 
function. The last term in Eq. (2.46) includes the effect of mass 
transfer. 


Eddy Viscosity 

The eddy viscosity model used is similar to that used by Tai [63], 
The expressions for the eddy viscosity are as follows: 


For attached flow, inner layer 



J 1 


u 

0.41 U 


8 = 1 + 0.0533 {e 


(l + 0.41 ^- + 0.5(0.41 -jj-) 2 )} (2.48) 


For attached flow, outer layer 


0.013 + 0.0038e 


-(6 /t )(dp/d0/15 


[' + 5.5 (-=) 6 ] 


U 6 Re 
e « 


(2.49) 


For separated flow, inner layer 


8=1+ 0.018U nRe fl - 

i e ® 1 V U ' J 


(2.50) 


For separated flow, outer layer 


0.01 3U 6 Re 
e ° 

[' * S.5(?f] 


(2.51 ) 


A substitution of — from Eq. (2.46) into the governing equations, 

T 

with U eliminated by evaluating U at the edge of the boundary 

p 

layer, results m three ordinary differential equations: 


A 11 

CN 

A 13 

A 21 

A 22 

A 23 

A 31 

A 32 

A 33 

B 3 

are as 

given i 



~(U ) 1 


B, 


T C 


1 


( 6 *) 

— 

B_ 


C 


2 


(U ) 


B, 


L e cj 


3 J 


(2.52) 


where A^ an 

The direct boundary-layer calculation corresponds to specifying 
U e and solving for 6* and U^. The inverse boundary-layer solution 


corresponds to specifying 6* and solving for U e and U^. Equation » 


i 
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(2.52) is reduced to the inverse form after rearranging the terms as 




(2.53) 


where 


C 


3 


are as given in the Appendix B. 


2.3.2 Initial Conditions 

The initial values for the viscous variables are evaluated, based 
on the Schlichting' s skin-f nction formula (65) for incompressible flow 
modified to include pressure gradient, as 



where m «= 4 and k = 0.0128. Starting with the mviscid edge 

velocity, both sides of Eq. (2.54) are evaluated until the required 
convergence condition is satisfied. For this value of U e , the 
frictional velocity is calculated using Eq. (2.55). 
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Chapter 3 

HETHOD OF SOLUTION 

The numerical procedures used to solve the direct and inverse 
boundary-layer equations are presented in this chapter. 

3.1 Direct Boundary-Layer Interaction 
In the direct approach the pressure distribution from the mviscid 
flow is specified and the boundary-layer quantities are evaluated. The 
mviscid flow over an airfoil is obtained by solving the potential flow 
equation in the entire flow field except m the shock region. In the 
shock region the Euler equations are solved using the information from 
the potential flow. This method is referred to as the hybrid method. 

In Sec. 3.1.1 the solution method to solve the potential flow 
equation and the Euler equations is discussed. A survey of the 

available viscous-inviscid interaction methods and the present method to 
achieve a continuous interaction between the mviscid and the viscous 
flow are presented in Sec. 3.1.2. Also, the description of the solution 
procedure for both flows is given. 

3.1.1 Solution of Inviscid Flow Equations 

The transonic potential flow equation, Eq. (2.3), subjected to the 
boundary condition, Eq. (2.5), is solved by the finite-difference 



J4 

relaxation scheme developed by Carlson [16]. In this method the 
governing equation is replaced by a nonconservative system of finite- 
difference equations and the system of equations are solved by a column 

t • 

relaxation technique. This procedure is adopted in this study because 
the difference equations are solved on a Cartesian grid. 


Solution of Euler Equations in the Shock Region 

It is necessary to solve the Euler equations m the shock region 
because of the rotational nature of the flow. It is also important to 
have a continuous interaction between the mviscid and viscous flows. 
In order to achieve this, the solution methods should be of the same 

type for both flows. For this reason the method of integral relations 

(MIR) is adopted to solve the Euler equations. Melnik and Ives [26], 

Holt and Mason [25], Sato [27], Tai [28], and Ram et al. [62] have used 

this method to solve the transonic inviscid flow equations for various 
flow conditions. Another advantage of using the MIR is its small 
computational requirement. 

The governing partial differential equations are reduced to 
ordinary differential equations by integrating Eq. (2.6) from the edge 
of the boundary layer to each strip boundary (Fig. 3.1) at some x 
location. In order to perform the integration, the integrand is 
approximated by an interpolation polynomial. For example, A in Eq. 
(2.6) can be approximated by 


N 

A = l a (x)(y 
k=0 * 




k 


(3.1 ) 
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Potential 

flow 

solution 

specified 


Line 3 (potential flow solution specified) 



Fig. 3.1 Integration scheme for the method of integral 


relation in the shock region. 



3b 


Using a second-order approximation for Eq. (2.6) the method can be 
implemented with three strips for desired accuracy. This process is 
( illustrated in Fig. 3.1. The integration domain is divided into two 
effective regions, which are denoted by strip boundaries (e, 1, 2) 

and (e, 2, 3). The base boundary e is set at the edge of the 
boundary layer. The flow conditions are specif iced by the potential 
flow solution on the uppermost boundary 3. First, the MIR is applied to 
determine the flow conditions along the boundary 2. Then, it is applied 
to the inner part of the flow field (e, 1, 2). 

The resulting ordinary differential equations for the mviscid 
external flow, reduced by means of the 2-2 strip integration scheme, 
associated with MIR, along the strip boundaries are 


dU 
e 

dx 

dV 
_ e 

dx 



F 
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e 
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G 
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(C 


1)(P t2 /P t2 ) ] ~ U 3 


< C - 1) ( p 2 





(3.2) 

(3.3) 

(3.4) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 
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P 

3 



(3.9) 


where 



and G are given m Ref. 41 . 


3.1.2 Shock/Boundary Layer Interaction and Viscous-Inviscid Coupli ng 

In the transonic flow regime the interaction between the boundary 
layer and the external flow is more important than in the subsonic or 
supersonic regimes. Also, the Reynolds number has a large effect on the 
aerodynamic characteristics as confirmed by Loving [66] through 
experiments m flight and in the wind tunnel. 

Experiments by Liepmann [29] and Ackeret et al. 130] indicated that 
the pressure rise m the boundary layer is much more gradual than in the 
external inviscid flow. When a normal shock impinges on the boundary 
layer, the disturbance propagates upstream through the subsonic portion 
of the boundary layer diffusing within a few boundary- layer thicknesses 
depending on the strength of the shock. If the pressure gradient is 
large enough, the flow may separate ahead of the shock. The pressure 
rise diverges the streamlines in the subsonic region generating 
compression waves in the supersonic region. In the case of the laminar 
boundary layer, the foot of the shock is thus smeared and a lambda shock 
appears. However, this does not necessarily happen in the case of the 
turbulent boundary layer because it can undergo a larger adverse 
pressure gradient than a laminar boundary layer. The displacement 
thickness increases considerably for a laminar boundary layer as 
compared to a turbulent boundary layer due to the shock. Also, the 
Reynolds number has a large effect on the laminar interaction but there 
is almost no effect on the turbulent boundary- layer interaction. 



After the shock the laminar boundary layer remains separated all 
along the airfoil due to the adverse pressure gradients encountered. 
Preferably, the coupling method should allow the downstream influence on 
the upstream. These types of viscous-inviscid coupling methods are 
termed strong interaction coupling methods and others are termed weak 
interaction coupling methods. 

Melnik [38] and LeBalleur [67] gave the recent state-of-art on the 
coupling of thin shear layer equations with the inviscid potential 
equations. In Bavitz's [68] method the effect of the wake is not 
included and an empirical correction near the trailing edge is used. 
Collyer and Lock [36] included the effect of a wake in their calcula- 
tions in the form of a normal velocity jump. However, they did not take 
the shock-boundary- layer interaction into consideration. Melnik et al. 
[38] have considered the trailmg-edge modelling but their method is for 
airfoils without any separation; also, the shock-boundary- layer 
interaction is not taken into account. 

Nandaman et al. [69] used Inger's [70] non-symptotic multi-deck 
analysis to predict a realistic pressure calculation in the shock 
region. They have used a solid displacement model with smoothing for 
interaction and the effect of a wake is not considered. Also their 
method is applicable to airfoils without separation. 

Wai and Yoshihara [71] considered an empirical model in the shock 
region to deal with separation. The interaction process is of semi- 
implicit nature. Updating the mesh periodically the curvature effect of 
the wake is taken into account by LeBalleur [67]. The interaction is 
acheived through a surface transpiration model. Klineberg and Lees [43] 
used the streamline angle at the edge of the boundary layer as a common 


i 



variable to allow for the continuous interaction between the boundary 


layer and the inviscid flow for supersonic external flow. Tai [41 ] , and 
Ram et al. [62] successfully employed this method for transonic flow. 
In the present study, this interaction model is used in the hybrid 
approach where both the viscous and inviscid solutions are obtained by 
using the integral methods. This facilitates a simultaneous interaction 
between the inner boundary layer and the outer inviscid flow. The 
common variable is given by the relation 

-1 V e i 

0 = sin ^ — - 6 (3.10) 

M a 
e,v e,v 

where V g ^ is the normal velocity component from the inviscid solution 

and M a is the magnitude of velocity from the viscous 

e,v e, v 3 

solution. The mass transfer between the two regions is allowed using 
the continuity equation 

* 

= tan 0 + ( 6 - 6 ) — ■ - in ( p U ) (3.11) 

ds ds e e 


where U g , the honzantal component of velocity at the boundary-layer 
edge is determined by the equation. 


U 

e 



(3.12) 


The strong interaction formulation can be applied to the attached, 
as well as separated laminar boundary layers, when it is applied to the 
attached flows, the boundary layer separates in a short distance. If 
the usual weak interaction formulation is applied, it is noted that the 
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separation occurs only at the shock. However, the strong interaction 
formulation involves another iteration process to determine the location 

of the shock influence point. In the forward portion of the airfoil, 

< * 

the weak interaction equations are sufficient to account for the 
interaction. 

Solution Procedure 

The transonic full potential equation is solved by the finite- 
difference scheme developed by Carlson [16], The shock location and the 
extent of the supersonic region is obtained from the Mach chart. This 
information is important to locate the shock influence point and to 
choose the strip boundaries to solve the Euler equations in the shock 
region. 

The initial displacement thickness and the velocity profile 
parameter are calculated using the procedure given in Sec. 2.2. In the 
forward portion of the airfoil the interaction between the boundary 
layer and the inviscid flow is considered to be weak. Therefore, the 
weak interaction formulation is applied to calculate the boundary-layer 
thickness and velocity profile parameter for a given pressure distribu- 
tion. The numerical integration of the boundary- layer equations is 
performed by a fourth-order Runge-Kutta method until the shock influence 
point is reached. 

Location of the Shock Influence Point 

To determine the shock influence point the strong interaction 
calculations are initiated at a number locations ahead of the shock. 


The shock location is determined from the mviscid solution 


The 
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displacement thickness (6*) and the velocity profile parameter 'a' 

are kept continuous when the switch is made from the weak interaction to 
strong interaction formulation. Also the velocity gradient at the edge 
of the boundary layer should be continuous to achieve convergence for 
the influence point. This is done by adjusting the streamline angle at 
the boundary-layer edge. 

In the inviscid flow region Euler equations are solved using the 
method of integral relations. The potential flow solution is taken as 
the initial condition along the vertical line and as the boundary 
condition along the outermost strip. The inviscid, as well as the 
viscous solutions, are obtained simultaneously. 

Usually the flow is separated shortly after the strong interaction 
equations are applied. If suction is not applied to keep the flow 
attached, the boundary- layer quantities are calculated based on the 
separated profiles. The integration continues downstream through the 
trailing edge and into the wake. At the trailing edge it is important 
to check the velocity gradient and adjust the streamline angle before 
continuing the calculation into the wake. 

The downstream boundary condition is satisfied for the correct 
shock influence point. The upper and lower surfaces are treated 
separately to compute the displacement thickness distributions. The 
velocity discontinuity at the trailing edge should be zero to satisfy 
the Kutta condition. This is checked by comparing the static pressures 
at the trailing edge from the upper and lower surfaces. The 
displacement thickness distribution is underrelaxed using a procedure 
discribed in Sec. 2.6. 





The airfoil is updated by adding the displacement thickness to the 
original airfoil coordinates and the inviscid flow is computed for the 
( updated airfoil with a new circulation accounting for the pressure 
difference near the trailing edge. With the new inviscid potential flow 
solution and the location of the shock, the procedure is repeated. This 
overall iteration process is continued until a specified convergence 
criteria on the displacement thickness is satisfied. 


3.2 Inverse Boundary- Layer Interaction 

In the inverse approach the boundary-layer thickness distribution 
is specified to avoid the separation point singularity. Semi-inverse 
coupling is used to couple the outer inviscid flow and the inner viscous 
flow. The solution method to solve the Euler equations m the entire 
flow field is presented in Sec. 3.2.1. The coupling method and the 
solution procedure for either a laminar boundary layer or a turbulent 
boundary layer are given in Sec. 3.2.2. 

3.2.1 Solution of Euler Equations Using Finite-Volume Approach 

The inviscid flow equations, Eqs. (2.29), are solved using the 
finite-volume approach developed by Jameson et al. [72). The 

discretization procedure decouples the spatial and time terms using the 
method of lines. The computational domain is divided into quadrilateral 
cells as shown in the Fig. 3.2, and a system of ordinary differential 
equations is obtained by applying Eq. (2.29) to each of these cells 
separately. This resulting system of equations is solved by the Runge- 


Kutta time stepping scheme 



Fig. 3.2 Computational grid Cor solving Euler 


equations by finite-volume approach. 
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The Runge-Kutta scheme has the advantage of allowing explicit time 
steps greater than a Courant number of one at the expense of evaluating 
additional functions at different stages. Whitfield et al. [73] has 
done some numerical experiments to find an optimum Courant number to 
perform the calculations. In this study, a four-stage Runge-Kutta 
scheme (with a Courant number of 2.8) is used; this is suggested in 
[74]. 

To suppress the oscillations near the shock and the stagnation 
points some external dissipation is added. The dissipative terms are a 
mixed blend of second-and fourth-order terms which are of third-order in 
smooth regions of flow and of first-order in the shock region. 

The convergence to steady state is accelarated by using a local 
time step (determined by the local Courant number) and by adding a 
forcing term that depends on the difference between the local and free 
stream values of total enthalpy. The fourth-order dissipative terms are 
needed to eliminate nonlinear instabilities when accelerating 
convergence using a local time step. 

3.2.2 Viscous-Inviscid Interaction Using Semi-Inverse Coupling 

When the Euler equations are used to compute the inviscid flow 
field the coupling requires momentum and enthalpy sources in addition to 
the mass sources. Johnston and Sockol [75], and Murman and Bussing [76] 
pointed out this information at about the same time. Thomas [74] 
followed the approach of Johnston and Sockol to achieve the viscous- 
inviscid interaction. Thomas modified the normal momentum relation by 
Rizzi [77] to include surface porosity for the pressure on the surface. 
In the present study the matching conditions are derived for the case 
with suction following Johnston and Sockol [75] . 
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Coupling Method 

The matching procedure adopted here is similar to that proposed by 
Johnston and Sockol 175] and is discussed here very briefly. The 
inviscid and the viscous solutions are matched on the surface. 

The Euler equations for steady two-dimensional flow can be written 


as 


where 


3F j)G 
3x + 3y 


0 


(3.13) 
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puv 
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u(e + p) 


v(e + p]_ 


The steady Navier-Stokes equations can be written as 


3f Jg_ 
3x + 3y 


0 . 


(3.15) 


In the defect formulation to be presented, the components of f and 
g are not needed. Integrating Eqs. (3.13) and (3.15) from y = 0 to 
y = 6, one obtains 

g 6‘ G 0 ‘ - * / F dy <3 -’ 6) 

- 5 o ■ - / £ ** (3 -' 7) 

By noting that, outside the boundary layer, the G and g vectors 
coincide, one can combine Eqs. (3.16) and (3.17) to obtain 
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3 6 

G o = g o + 17 / (f - F)dy {3 * 18) 

0 

, Now representing the Navier-Stokes solution by a composite function of 
the type 

f=f = F + f - F (3.19) 

c b 0 


where f^ is the boundary-layer solution, Eq. (3.18) reduces to 


Vo 


JL 

3x 


6 

/ (p 0 - f b )d » 


(3.20) 


It should be noted that only the value of F at the wall is needed and 
specific variation of F in the boundary layer is not necessary. 

From Eq. (3.14) we can evaluate the values of the vector G. The 
first term of G is expressed as 


(pv) o = Vo + &T / t (pu) o - Pb u J dy 


(3.21 ) 


Using the definition of displacement thickness Eq. (3.21) can be written 
as 


<P V) 0 


d 

dx 


[ ( pu ) ^ 6* ] + (P b v b ) 0 


(3.22) 


The second term of the vector G is written as 

g 6 2 2 

(puv) - v Vo - T 0 + / [(pu + P) 0 " °b + P b)l dy 

0 

By considering the pressure from the boundary- layer solution to be equal 
to the pressure from Euler solution and using the definition for the 
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momentum thickness, Eq. (3.22) can be written as 

( Puv) q = “ T 0 + ^ [ ( P u) 0 ( 6 + 6) ] (3.23) 

To evaluate the third term of G the pressure at the wall was 
determined by Thomas [74] using the Rizzi's normal momentum relation 
[77]. The surface porosity term was included in this analysis. 
Different interpolation relations to obtain pressure on the body from 
the adjacent cells are suggested by Thomas [74]. In this study some of 
these relations are used to check the convergence rate and the accuracy 
of the results. The fourth term of vector G can be written as 


({e + P)v) Q = [(e + P)v - uT] fa 


+ 'b ! U {e + P )u ] 0 " l (e + p )u ] b l dy 

0 


(3.24) 


With the approximation that the total enthalpy from the boundary layer 
is equal to the Euler solution value, Eq. (3.24) becomes 


[(e + p)v] Q = [(e + p)v] b + H q — [( pu) Q 6 ] 


pv Ho = (pH v) b+ H 0 — [(pu) 0 6 ] 


(3.25) 


This is an identity. As can be noted from (3.22) and (3.23), only the 
wall values are needed from the Euler solution and they can be obtained 
easily once the boundary- layer solution is obtained. 
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Wake Relation 

In the wake the boundary-layer quantities are evaluated along the 
,wake center line for upper and lower surfaces separately. The 

difference in the mass, momentum, and energy fluxes is applied as the 
boundary condition to the inviscid flow. In this study, neither the 

wake curvature effect nor the strong interaction near the trailing edge 
are taken into consideration. The computational mesh needs to be 
recomputed periodically to take the curvature effect into account. 


Viscous-Inviscid Coupling 

The viscous-mviscid interaction is achieved through a semi-inverse 
coupling. This technique was developed by Carter [78] for subsonic 
flows and has been used for transonic flows by Whitfield et al. [73], Le 
Balleur [67], and Thomas [74]. The inviscid algorithm is advanced 100 
time steps to obtain an approximate pressure distribution around the 
airfoil. Then the inverse boundary- layer equations are solved using a 
specified displacement thickness distribution. Initially, this 
distribution is that of a flat plate. The velocity at the edge of the 
boundary layer is calculated at all the grid points on the airfoil and 
along the wake center line. Then the semi-inverse coupling compares the 
velocities at the edge of the boundary layer obtained from viscous and 
inviscid solutions. Then the initial distribution of the bound ary- layer 
thickness is updated using the relation 


(« ) 



♦ «(£*- 1 )] 

e, i 


new 


(3.26) 
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where u is the boundary-layer edge velocity computed from the 

e i v 

inverse intergral method, u . is the edge velocity from the inviscid 
solution, and u is a relaxation parameter. 

By using the new boundary-layer thickness distribution, the 
boundary condition to the inviscid flow is computed by Eqs. (3.22) - 
(3.25). Then, the inviscid algorithm is advanced in time with this new 
boundary condition near the wall. The inverse boundary-layer equations 
are solved after every 20 inviscid cycles from then on and the boundary 
conditions are updated. When there is no appreciable change in the lift 
or the boundary- layer thickness distribution, the solution is considered 
to be converged and the calculations are ended. 

Suction 

The reduction in drag is an order of magnitude from the turbulent 
boundary layer to laminar boundary layer. Therefore, large extents of 
laminar flow are desirable to increase the aircraft performance. 

The laminar boundary layer cannot undergo large adverse pressure 
gradients and the flow separates resulting in rapid thickening of the 
boundary layer j this increases the drag. To avoid this undesirable 
effect, suction can be applied before the flow separation point to keep 
the flow attached all along the airfoil. 

Location of Suction 

The velocity profile parameter 'a' is proportional to the 
velocity gradient near the wall. The value of ’a' is zero at 

separation. A specified amount of suction is applied when the value 
of 'a' falls below certain value depending on the airfoil under 
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consideration. The value of the suction parameter is increased at the 
next location if the value of 'a' is still below the specified 
, value. When there is a negative pressure gradient or when the value 
of 'a' is higher than a specified value, the suction parameter is 
decreased. 

With the above procedure, the flow separation is avoided and 
attached flow condition is maintained all along the airfoil. The amount 
of suction thus calculated is close to the minimum amount. 
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Chapter 4 


RESULTS AND DISCUSSION 


Results have been obtained for different airfoils by employing the 
direct and inverse boundary- layer procedures for laminar and turbulent 
flows. Specific results were obtained for various cases with and with- 
out suction in the boundary layer. Results of calculations for the 6% 
circular arc and LFC-73-06-135 airfoils are presented in Sec. 4.1. 
These results were obtained for laminar flows using the direct boundary- 
layer equations. In Sec. 4.2 results are presented for the King Cobra 
airfoil and the modified NACA 66-012 airfoil. The flow conditions are 
selected such that the results can be compared with the experimental 
data. Also results for DESB-154 and LFC-73-06-135 airfoils are 
presented for attached flow conditions. Results for the laminar inverse 
boundary- layer equations coupled with the Euler equations are presented 
in Sec. 4.3 for the NACA-0012 airfoil and the RAE-2822 airfoil. Using 
the same approach, results were calculated for the turbulent boundary- 
layer flows with and without suction; these are presented in Sec. 4.4. 

4.1 Laminar Direct Boundary- Layer Solutions Without Suction 

Results of calculations at supercritical freestream Mach numbers 
are presented for a 6% circular arc and LFC-73-06-135 airfoils. Flow 
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conditions were chosen to enable comparisons with available experimental 
data . 

The viscous results were calculated in terms of boundary-layer 
quantities in a transformed incompressible plane. Figure 4.1 gives the 
boundary-layer displacement thickness throughout a 6% circular arc 

4 

airfoil at M - 0.868 and Re = 6.9 x 10 , which agrees very well 

co ao 

with similar results reported elsewhere [79] . The thickening of the 
boundary layer in the forward portion follows a similar trend as that 

found by Schubauer, using the Karman-Polhausen method [80]; however, the 

* 

present method gives a far more realistic 5 distribution pattern in 
the rear portion. 

The corresponding distribution of the boundary- layer shape factor 
H and the velocity profile parameter, a, are presented in Figs. 4.2 and 
4.3, respectively. The boundary layer is practically, but not exactly, 
of Blasius type in the leading-edge region and varies slightly through- 
out the forward portion of the airfoil. It remains unseparated through 
the embedded supersonic region although the viscous-inviscid interaction 
becomes strong after x * 0.395. The separation point is found when 
a = 0 which corresponds to zero shear stress at the wall. 

The boundary layer remains separated over the rear of the airfoil 
where small adverse pressure gradients are generated by continuous 
compression of the outer subsonic flow. This is a physical feature of 
the transonic viscous-inviscid interaction since by compression the flow 
ought to return almost to the free-streara value downstream. After the 
trailing edge there is a wake reversed flow and then a forward flow to 
match the downstream conditions. The location of the rear stagnation 
point agrees well with that found by Klineberg and Steger [39], and by 
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Fig. 4.3 Velocity profile parameter for a 6% circular arc 

4 

airfoil at M ** 0.868 and Re « 6.9 x 10 . 
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Tai [79] under similar flow conditions. Computed results for the 6% 
circular arc airfoil compare very well with the laminar experimental 
data of Collins and Krupp [81] as presented m Fig. 4.4, not only in 
pressure distribution but also in separation point. The small 

difference in free-stream Mach number and Reynolds number between theory 
and experiment was selected deliberately to offset wind-tunnel inter- 
ference effects [82]. 

There is no referenceable experimental data available at the 
present time for supercritical airfoils for which the boundary layer 
remains laminar over most of the airfoil. Of course, plenty of 
experimental data are available for many airfoils with turbulent 
boundary layers. For this reason it was decided to compare the theo- 
retical predictions of the present method to that of some other existing 
methods, such as Carlson's TRANSEP [83], to judge its reliability and 
accuracy. Based on the experience with the ongoing swept supercritical 
LFC airfoil experiment m the Langley 8-Foot Transonic Pressure Tunnel 
[84], it is expected that the flow will remain laminar over an extensive 
chordwise length of the LFC-73-06-1 35 airfoil for M = 0.75, Re = 

00 Q 

6 o 

8x10, and a = -0.09 . Therefore, this airfoil and these flow 
conditions were chosen to calculate viscous results using both TRANSEP 
and the present method, along with the inviscid results using Carlson's 
TRANDES [85]. Figures 4.5 and 4.6 show the pressure distribution 
obtained from these three methods. The viscous results obtained from 
these two methods are in very good agreement with each other except for 
a few deviations which were expected. For example, the c^ values for 
the upper surface obtained from TRANSEP are slightly lower (i.e., more 
negative) than those from the present method in the region lying between 
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Fig. 4.5 Pressure distribution on the upper surface of a swept 

supercritical LFC-73-06-1 35 airfoil at = 
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6 Pressure distribution on the lower surface of a swept 

supercritical LFC-73-06-1 35 airfoil at = 0.750, 
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39% through 76% of the chordwise length. The starting position of the 
adverse pressure gradient is located at 50% and 46% of the chord, 

respectively, as predicted by TRANSEP and the present method. The 

present method predicts flow separation at x = 0.665, whereas the flow 
remains attached over the airfoil according to TRANSEP. Because of the 
strong interaction formulation modeled in the present method, the flow 

separates, whereas it remains attached due to weak interaction formula- 
tion modeled in TRANSEP. Also, the displacement thickness of the 

laminar boundary layer increases more rapidly through the shock than 

that of turbulent boundary layer [86] and this may cause laminar 

separation. The pressure rises more rapidly for turbulent than for 
laminar boundary layers [87,88]. From 80% of the chord to the trailing 
edge, the values predicted by TRANSEP are higher (i.e., more 

positive) than those obtained from the present method. This is true 
because the boundary layer over the last 20% of the chord is definitely 
turbulent, whereas m the present method the boundary layer over the 

entire chord length always remains laminar. Thus, the theoretical 
results calculated from this method conform very well to the findings of 
[ 86 - 88 ] . 


4.2 Laminar Direct Boundary-Layer Solutions With Suction 

Results of calculations are presented for modified NACA 66-012 
[89], DESB 154 [90], King Cobra [91], and LFC-73-06-1 35 [92] airfoils. 
The airfoil NACA 66-012 was chosen to enable comparison of the calcu- 
lated result with experimental data available from [89] . The DESB 154 
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and King Cobra airfoils retain laminar flow over an extensive chordwise 
length (approximately 70% and 65%, respectively). These two airfoils 
were chosen because of extensive laminar flow which will be more 
appropriate to test the accuracy and reliability of the present method. 
Furthermore, experimental data available for a King Cobra airfoil for 
the no suction case permits comparison of the computed results under 
similar conditions. Since the method has been developed for super- 
critical airfoils, LFC-73-06-1 35 airfoil designed at NASA Langley was 
selected for comparison. 

The viscous results were calculated m terms of boundary- layer 
quantities in a transformed incompressible plane. Figure 4.7 gives the 
boundary- layer displacement thickness on the upper surface of a natural 
laminar flow airfoil (DESB 154 at M = 0.4, Re = 10 x io 6 , and a = 

ao 

O 

-0.97 ) for different values of the suction parameter e. It was found 
that e = 0.015 was the minimum amount of suction that kept the flow 
attached all the way to the trailing edge. The suction was started at 
65% of the chord length and was maintained up to the trailing edge. It 
should be pointed out here that the flow separates at about 70% of the 
chord length for the same flow conditions if the suction is not applied, 
i.e., e = 0, as shown by the present theoretical computations as well as 
by those in [90] . It is further seen that the thickening of the 
boundary layer and hence, the overall viscous effects can be controlled 
easily by varying the magnitude of the suction parameter e. 

The corresponding distribution of the velocity profile parameter 
a is presented in Fig. 4.8. There are minimal fluctuations in its 
values until the point of separation (in the absence of suction) is 
approached when it suddenly increases and remains high or drops down 






Fig. 4.8 Velocity profile parameter for DESB 154 airfoil at 

M = 0.4, Re = 10 x 10 6 , a = -0.97° for different value 
00 

of the suction parameter e. 



64 


depending on the amount of suction applied. The lowest value, a = 0, 
occurs at the point of separation and corresponds to zero shear stress 
at the Wall. 

Figure 4.9 shows the pressure distribution over a DESB 154 airfoil 
obtained by using the Carlson's TRANDES (Inviscid) method and the 

present method with suction. These theoretical results compare very 
well as expected. The values for the upper surface obtained from 

the TRANDES (Inviscid), when there is no boundary layer, are slightly 
lower (i.e., more negative) than those from the present method with 
suction, when there is a very thin boundary layer. The point of sudden 
rise (i.e., more positive) in value occurs at 0.71c and 0.68c, 

respectively. 

Figures 4.10 and 4.11 give the pressure distribution over the King 
Cobra airfoil without and with suction, respectively. The computed 

results in Fig. 4.10 compare very well with the experimental data 
reported in (91], not only in pressure distribution but also in 

separation point. The suction for the King Cobra airfoil was started at 
60% of the chord length and was maintained up to the trailing edge. It 
should be mentioned here that the flow separates at about 65% of the 
chord length for the same flow conditions if the suction is not applied, 
i.e., e = 0, as shown by the theoretical computations as well as 

experimental data. 

Figure 4.12 shows the pressure distribution over a swept LFC 
airfoil (modified NACA 66-012) obtained from the present method with 
suction. The computed results compare very well with the experimental 
data, also with suction, reported in (89], 
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Fig. 4.9 Pressure distribution on the upper and lower surfaces ot 
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DESB 154 airfoil at M = 0.4, Re = 10 x 10 , o = -0.97 , 
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E ** 0.015. 
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4.10 Pressure distribution over King Cobra airfoil at 
M « 0.4, a = 1°, Re = 10 x io 6 and e » 0 . 
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Fig. 4.11 Pressure distribution on the upper and lower surfaces of 

6 o 

King Cobra airfoil at M = 0.4, Re = 10 x io , a = 1 , 
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C = 0.199, e= 0.002 for upper surface obtained from the 


present program. 
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Fig. 4.12. Pressure distribution on the upper surface of a modified 
NACA 66-012 airfoil at M = .12, o *= 0°, Re = 14 x IQ 6 , 


and .003 4 z < .015. 
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There does not exist any computational method which solves strong 
mviscid-lammar viscous interactions by solving Euler equations in the 
shock-dominated region and where the boundary layer is optimally sucked 
such that boundary-layer instability and separation do not occur. The 
pressure distributions over the LFC-73-06-1 35 supercritical airfoil 
calculated from Carlson's TRANDES (Inviscid) and the present method are 
compared to show the effect of suction and the formulation of the 
viscous-inviscid model. The flow conditions were chosen to assure the 
existence of a shock on the upper surface of the airfoil. In the 

inviscid analysis, Fig. 4.13, the shock appears at 77% of the chord 
length. In the mviscid-viscous analysis with suction, Fig. 4.14, the 
shock becomes much weaker, moves upstream to 0.63c and the boundary 
layer remains attached. Figure 4.15 shows that C p values calculated 
from the present method are consistently higher (l.e., more positive) by 
as much as up to 24.3%. This is expected due to effect of viscosity and 
the change of entropy across the shock which are accounted for in the 
present method. The results mentioned here for LFC-73-06-1 35 airfoil are 
in agreement with the findings of the swept supercritical LFC airfoil 
experiments conducted in the Langley 8-Foot Transonic Pressure Tunnel 
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Fig. 4.13 Pressure distribution over a swept supercritical 
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LFC-73-06-1 35 airfoil at m = 0.780, a = 0 , Re = 10 x 10 , 

and e = 0 obtained from Carlson's TRANDES (Inviscid). 
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Fig. 4.14 Pressure distribution over a swept supercritical 
LFC-73-06-1 35 airfoil at M = 0.780, a = 0°, 
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Re = 10 x io 6 , e = 0.035 on the upper surface 
obtained from the present program. 
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Fig. 4.15 Pressure distribution on the upper surface of a 
swept supercritical LFC-73-06-1 35 airfoil at 
= 0.780, a= 0°, and Re = 10 x io 6 . 
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4.3 Inverse Laminar Boundary-Layer Solutions 

Results are presented in Fig. 4.16 for an RAE-2822 airfoil at 
subcntical condition. The free stream Mach number is 0.6 and the angle 
of attack is 2.57°. The inviscid flow calculations are performed by 
solving the unsteady Euler equations [72]. The viscous-inviscid 
interaction was started very close to the leading edge. A small amount 
of suction was required to keep the flow from separating near the 
leading edge. The pressure peak near the nose causes early transition 
to turbulent flow if suction is not applied. After the pressure 
minimum, the flow is continuously decelerated to match the downstream 
flow of the airfoil. The laminar flow is stable in this region and the 
flow remains attached all the way to the trailing edge. There is no 
need to apply suction after the nose peak. The boundary layer is 
thicker than in the case of a turbulent flow and the values, are 

less negative along the airfoil. 

Figure 4.17 illustrates results for the supercritical conditions. 
The free stream Mach number is 0.73 and the corrected angle of attack is 
2.78°. This case corresponds to case 9 of the experiments conducted by 
Cook et al. 194]. In the experiments, the flow was tripped very close 
to the leading edge to produce a turbulent boundary layer. For laminar 
boundary layers the flow would separate at about 45% of the chord 
length. The suction was applied before the separation point to keep the 
flow attached. The amount of suction applied depends upon the velocity 
profile parameter a. The value of the suction parameter is increased 
or decreased based on whether the value of a is below or above a 
specified limit. This was an effort to apply only the minimum suction 
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Pressure distribution on the upper surface of RAE-2822 
supercritical airfoil at = 0.6, a = 2.57°, and 
Re «* 6.5 x 10 6 . 


Fig. 4.16 


75 



Fig. 4.17 Pressure distribution on the upper surface of RAE-2822 
supercritical airfoil at M = 0.73, a = 2.78°, and 

OO 

6 

Re = 6.5 x 10 , 
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required rather than an arbitrary amount. The strength of the shock is 
reduced and moved forward because of the laminar boundary layer as 
•compared to the mviscid flow calculations. 

For NACA-0012 airfoil extensive experiments were conducted by 
Harris for subcntical to supercritical Mach numbers and for different 
values of the angle of attack [95] . However, in all these experiments 
the boundary layer has been tripped close to the nose to produce a 
turbulent boundary layer flow. Present results for transonic laminar 
flow are compared with these experimental data for a qualitative com- 
parison. Results presented in Fig. 4.18 are for a Mach number of 0.758 
and an angle of attack of 3.06°. The laminar viscous-inviscid inter- 
action was started just ahead of the stagnation point at about 1% of the 
chord length on both the upper and lower surfaces. The location of the 
shock foot is at about 52% of the chord length according to the mviscid 
flow results. From the viscous-inviscid interaction calculations, it is 
noted at about 48% of the chord length. The suction was applied at 
47.48% to maintain the laminar attached flow. A large amount of suction 
was required up to about 55% of the chord and after the shock, a rela- 
tively small amount of suction was sufficient to keep the flow attached. 
In the experiments the turbulent flow was separated at 34% of the chord, 
the shock strength was reduced considerably and moved forward in com- 
parison with the present results. The suction was applied to produce a 
turbulent, attached boundary layer flow over this airfoil under similar 
conditions and those results are presented in the next section. The C 

P 

values of both the experimental investigation and the present method are 
close to the mviscid results in the rear part of the airfoil where the 
flow has to match the subcritical conditions in the downstream region. 
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Fig. 4.18 


Pressure distribution on the upper surface of NACA-0012 
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airfoil at M = 0.758, a = 3.06 , Re = 3.0 x 10 . 
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Results are presented in Fig. 4.19 for an NACA-0012 airfoil at high 
angle of attack conditions. The free stream Mach number is 0.601 and 
Ithe Reynolds number based on the chord length is three million. A 
correction to the angle of attack was applied to compare the results 
with the experiments [95]. It is evident from the Fig. 4.19 that there 
is a pressure peak in the forward portion of the airfoil and there is a 
pressure rise at about 17-20% of the chord. The experimental values are 
closer to the inviscid flow results in the peak region and the pressure 
rise region except near the end of the shock. The viscous-mviscid 
interaction was started at about 1% of the chord. Suction was applied 
to maintain laminar attached flow. After the shock not much suction is 
required. In comparison to the inviscid flow, the shock strength is 

reduced and the pressure wiggle is reproduced at the end of the shock. 
The coefficient of lift is 0.761 as compared to the turbulent value of 
0.847. 


4.4 Turbulent Boundary-Layer Results With Suction 

Results are presented in Fig. 4.20 for the RAE-2822 airfoil at 
subcritical conditions. The free stream Mach number is 0.6 and the 
angle of attack is 2.57°. A correction to the angle of attack is made 
following the suggestions of Cook et al. [94]. The results are compared 
with the experimental data of Cook et al. [94] and with the theoretical 
results of Thomas [74] . The Reynolds number based on the chord length 
is 6.5 million. The viscous-inviscid interaction was started at about 


15% of the chord on the lower surface and at 18% of the chord on the 
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upper surface. In the experiment the boundary layer has been tripped at 
these locations. The initial displacement thickness values correspond 
to the experimental data. The calculations were performed on a C-type 
mesh using 128 points along the airfoil and 30 points away from the 
airfoil. The boundary layer was iterated for every 20 cycles of the 
inviscid calculations (after the first 100 cycles). In the forward 
portion of the airfoil, where there is a pressure peak, results obtained 
from the present calculations agree closely with the experimental 
values. The turbulent flow remains attached until the trailing edge and 
there is no suction applied in this case. The lift coefficient, C^, 
values were found to be 0.71 and 0.65 for inviscid and turbulent flows, 
respectively. The results for the RAE-2822 airfoil for supercritical 
conditions are illustrated in Fig. 4.21. The free stream Mach number is 
0.73, the corrected angle of attack is 2.78°, and the Reynolds number 
based on the chord length is 6.5 million. These conditions correspond 
to case 9 of Ref. 94. These calculations were also obtained using a C- 
typ>e mesh with 128 x 30 points. The grid was highly stretched away 
from the airfoil until the change in the lift coefficient was small 
[74]. The inviscid lift coefficient showed no appreciable change after 
600 time steps. The extent of reduction in the maximum residual is of 
fourth-order. After the first 100 time steps the boundary layer was 
interacted for every five time steps. Although this increases the 
computational time, the results were more accurate. Because of frequent 
updating, the viscous-inviscid interaction process is closer to the 
strong interaction conditions. The interaction was started at about 18% 
of the chord with the initial displacement thickness values from the 
experiment [94]. For comparison, the equilibrium dissipation model of 
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Thomas [74] was used. The pressure coefficient values are compared in 
Fig. 4.21. In the forward portion of the airfoil the pressure 

coefficient values are slightly lower with the present model. In this 
area results using the equilibrium dissipation model of Thomas [74] are 
closer to the experimental data. The agreement between the results of 
the two methods is very good except close to the shock foot. The 
values are underpredicted near the shock and in the rearward portion of 
the airfoil. The lift coefficient values are found to be 0.97,0.90 and 
0.91 for the inviscid, equilibrium dissipation, and present methods, 
respectively. 

Results are presented for a supercritical case with separation at 
about 60% of the chord m Fig. 4.22. This case corresponds to case 10 
of Ref. [94]. The free stream Mach number is 0.75 and the Reynolds 
number of the flow is 6.2 million. The corrected angle of attack is 
2.81°. The viscous-inviscid interaction calculations were started at 
about 18% of the chord on the upper surface and at 15.5% on the lower 
surface. The boundary layer has been tripped to become turbulent at 
these points m the experiment. The present results are compared with 
the experimental data of Cook et al. [94]. The pressure distribution m 
Fig. 4.22 indicates that the agreement with the experiment is good in 
the acceleration zone where the boundary layer is thin. In the 
experiment the flow separation was observed near the foot of the shock 
between 62-72% of the chord length. In the present method suction was 
applied near the shock foot to keep turbulent attached flow conditions 
on the airfoil. The displacement thickness values are compared in Fig. 
4.23. Because of suction, the displacement thickness does not increase 
as rapidly as the experimental values in the pressure rise area. In the 
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o Turbulent flow (present) 



Fig. 4.23 Displacement thickness distribution on the upper surface of 
RAE-2822 supercritical airfoil at M = 0.75 a= 2.81° 

ao r 

, 6 
and Re =6.2 x 10 . 
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experiment the wake is much thicker due to turbulent separation compared 
to the turbulent attached flow results by the present method. Because 
of this reason the values are higher near the trailing edge. 

The upper surface pressure distribution is presented in Fig. 4.24 
for the NACA-0012 airfoil at subcntical conditions. The free stream 
Mach number is 0.601 and the corrected angle of attack based on the 
linear theory is 3.19° 195]. The Reynolds number of the flow is three 
million. The results are compared with the experimental data of Harris 
[95] . In the experiments the flow was tripped at about 5% of the chord 
length. The boundary-layer interaction was started corresponding to the 
experimental data. The numerical results obtained in this study and by 
Thomas [74] agree closely with the experimental data. However, better 
agreement is noted between the experiment and the present method in the 
forward portion of the airfoil. In the rearward portion, the pressure 
distribution is very close to the inviscid case. 

The pressure coefficient results are presented in Fig. 4.25 for a 
free stream Mach number of 0.758 and for an angle of attack of 3.06°. 
Suction was applied at about 30% of the chord to keep the flow from 
being separated. In the experiment there was no attempt to apply any 
suction. The flow might have separated at about 35% of the chord and 
reattached after the shock. For the above conditions, the computer 
program using the equilibrium dissipation model failed to produce any 
results because of the extensive separation of the flow. Also the 
present method does not work if the suction is not applied before the 
separation occurs. The agreement between experiment and the present 
method is good up to about 30% of the chord. The turbulent boundary 
layer thickens rapidly after that point. The flow is separated and the 
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Fig. 4.25 Pressure distribution on the upper surface of NACA-0012 
airfoil at m = 0.758, a * 3.06°, and Re = 3.0 x 10 6 . 
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shock front is moved towards the leading edge and also, as observed m 
the case of laminar flow, the strength of the shock is reduced. By 
applying the suction before the separation point, the flow separation 
was delayed by about 20% of the chord. The pressure rise is much more 
gradual that it would be in the case without the suction. After the 
shock the computational results are close to the experimental data. 
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Chapter 5 


CONCLUSIONS 


The viscous/inviscid interaction over transonic airfoils with and 
without suction is studied. Two approaches are considered to achieve 
the coupling between the viscous and inviscid flows. The first approach 
is a direct approach and is referred to as the hybrid method. In the 
second approach the entire inviscid flow field is investigated by 
solving the Euler equations using finite volume technique. The viscous 
flow is coupled to the inviscid flow using surface transpiration 
condition. 

The interaction process m the hybrid method is continuous, and 
since all the dependent variables are calculated simultaneously, the 
convergence is faster and the solutions are more accurate. Using this 

method, flow over a 6% thick circular arc airfoil at M = 0.868 and 

o 

a = 0 is studied. In the forward portion of the airfoil it was 
sufficient to apply weak interaction formulation and the strong inter- 
action equations are applied near and downstream of the shock. The 
separation was predicted at 70% of the chord and is in complete agree- 
ment with the experimental data of Collins and Krupp [81 J. The pressure 
distribution calculations using laminar separated velocity profiles show 
a good agreement with the experimental values [81]. 
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To avoid the Goldstein singularity near the separation point, the 
boundary- layer equations and the initial conditions are derived in an 
inverse ' form to obtain regular solutions for the flows with small 
separations. It is important to apply minimum amount of suction that is 
required to have attached flow on the airfoil. A method is developed to 
achieve this by varying the suction parameter based on the velocity 
profile parameter value. The suction distribution thus obtained is 
considered to be close to the optimum value. At subcntical conditions 
the present solutions are compared with the experimental data [89-91) 
and the agreement is excellent for NACA 66-012, DESB 154 and KING COBRA 
airfoils. These comparisons for subcntical airfoils establish 
confidence in the suction velocity profiles that are obtained using the 
small perturbation theory and the state variable approach of Forbrich 
[60]. 

Results are obtained for RAE-2822 and NACA-0012 airfoils at super- 
critical conditions. These results indicate that a small amount of 
suction is required to avoid flow separation near the pressure peak at 
the leading edge. Most of the suction requirement is to maintain 
attached flow conditions in the rear part of the airfoil or near the 
shock. Application of larger amount of suction than required was found 
to have destabilizing effect on the boundary layer. Laminar* flow 
separation reduced the shock strength considerably and shock is moved 
forward in comparison to the laminar attached flow with suction. 

The results from the present study are in good agreement with the 
theoretical as well as the experimental data for attached flow condi- 
tions with turbulent boundary layer. At supercritical conditions the 
boundary layer tends to separate and application of suction has been 
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considered to maintain attached flow. Although the eddy viscosity model 
used is valid for both attached and separated flow conditions, the flow 
separation calculations could not be performed satisfactorily, perhaps 
due to the presence of suction. If the experimental data is available 
for turbulent separated flow with suction, the present method could be 
extended to separated flow conditions. For laminar as well as turbulent 
boundary layers the displacement thickness is small compared to the 
corresponding separated flow conditions and the wake thickness is much 
smaller. The strong interaction near the trailing edge has to be 
considered if the flow separates near the shock. 

The boundary-layer integral method coupled with the method of 
integral relations gives a computationally inexpensive solution for 
transonic laminar viscous-inviscid interaction over airfoils. However, 
this method requires man-machine interaction and the solution can not be 
obtained in one computer run. The inverse boundary-layer approach 
obtains the flow solution as well as the suction distribution to keep 
the attached flow on the airfoil in one run. However, the boundary 
layer equations are not of strong interaction type. Since the shock 
influence point has to be determined in an iterative process the 
computational requirement is very high. The vectorized version of the 
Euler solver for the inviscid flow could be used in a further study to 
include the shock influence point iteration so that the strong interac- 
tion model is incorporated. 
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Appendix A 


LAMINAR BOUNDARY-LAYER EQUATIONS WITH SUCTION 


The laminar boundary-layer equations (2.15) can be expressed 
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Appendix B 


TURBULENT BOUNDARY-LAYER EQUATIONS WITH SUCTION 


The turbulent boundary-layer equations are given by Equation (2.52) 


as 


1 1 

A 12 

A 13 

‘21 

A 22 

A 23 

31 

A 32 

A 33 


(U ) „ 


B, 

T C 


1 

(6*) 

_ 

B_ 

C 


2 

(U ) 


B 

e C 


3 _ 


(B.1 ) 


where 


6 

A = J FIRST dn 

0 

6 

A = / SECOND dn 

2 0 

6 

A = j THIRD dn 

13 0 

6 

A 21 = / (FIRST) FACT dn 

6 

A - J (SECOND) FACT dn 

22 0 



105 


FIRST » 


6 

a„, - / (third) fact dn - 0.5 u 
23 0 

6 

A = J (FOURTH) dn 

0 

6 

A = / (FIFTH) dn 

0 

6 

A = / (SIXTH) dn 
0 


FACT = 


n 6 

w - v e - ! dn + / dn 
0 e 0 e 


B , ■ 


v °’ e“) 


Re 6U '•9n'w 

00 0 


B 2 ‘ 


6 ^ d " 


B = V - V 
3 w e 


d6 

C 1 “ B 1 A 12 dC 


C = B - 
2 3 


d6 

A__ ~ — 
32 d? 


U ( —0 .5(1 — C2 ) J [2.5 log (C6) + 5. I) 2 V 

T 1 W 


_ pC 

+ 2.5 log(C6) - (3.39 DELP + 5.1 )e + 5.1 ] 


+ [(2.5 log (C7) + 5.1) V w + 2.5 log (C7) 


- (3.39 ETAP + 5.1 )e + 5.1 ] 



1 06 


j 


+ 0.5(1 - C2)U 1 |+ 0 . 5 (C2 - UU.J- 


5. 0(2. 5 log (C6) + 5 . 1 ) 
C6 


C9 


4- 0.37(3.39 DELP + 5.1)e _C5 C9 6 - 3.39 e' C5 C9 6 


2.5 C9 6 
C5 


5.0 n C9 V 


]+ U T [ (2.5 log (C7) + 5.1) 


+ 0.37(3.39 ETAP + 5 . 1 ) n C9 e _C5 - 3.39 n e _C5 C9 


2 c n pg 2 

+ ] - 0.5(1 - C2 ) [(2.5 log C6 + 5.1) 


-C5 

+ 2.5 log C6 - (3.39 DELP + 5.1 )e +5.1) 


+ (2.5 log C7 + 5.1) V w + 2.5 log (C7) 


- (3.39 ETAP + 5.1 )e _C5 + 


5.1]J 


w 


j-0.5 Cl it [ ( 2 . 5 log C6 + 5. 1) 2 

pc 

+ 2.5 log C6 - (3.39 DELP + 5.1 }e + 5.1 ]/ 6 


5. 0(2. 5 log C7 + 5.1 ) |u I V 

1 T W 

C7 v 


0.37(3.39 + 5.1 - e~ C5 luj 


2.5 (U t | 0.5 Cl if ul 


C7 v 


5. 0(2. 5 log C6 + 5.1) C9 V 6 
0.5 U t (C 1 6 - Tin) [' gjj — ] 


V 6 
w 



107 


+ 0.37(3.39 DELP +5.1 - ) e ~ C5 C9 6 

U • J r 

2.5 C9 , 

r 2 

+ 0.5 [(2 .5 log C 6 + 5.1) V w 

—C5 

+ 2.5 log C 6 - (3.39 DELP + 5.1 )e + 5.1 ] 

— C5 C 8 2 

• [ci 6 - nn]/it + 7.305 e [(e (0.8556 n |u 

+ 0.8556 vU 2 t )log 2 C7 + e C 8 [l.7797 U 2 (n |u | + v)log C7 
+ 1 .7811 r*e C8 U 2 |u I ]v 

T 1 T* J W 

+ e C8 [0.3423 (J 2 ( nU^ + v) ]log C7 

3 C 8 2 2 

+ 1 .2543 n |U t | + e (0.3560 nU^ |u - 5.277 vU^) 

+ 5.277 VU 2 ]/(u 2 |U x | ) 

-C 8 T C 8 2 2 , , 2 

- 7.305 e (0.8556 e v U T l u T l lo 9 c7 

+ e C8 (1 .7797 v 2 U 2 |U t | - 1.71125 nv U 4 )log C7 

- 1 .7797 ne 08 v U 4 V 

T j W 

v|(J | 

+ 0.3423 e C8 v 2 U 2 |u |log( 5 -) + l) 

T T nu 



108 


SECOND 


C8 2 2 

- 0.3423 e v U T i u T ! l°g( v juj) 


PR p p p p 

+ e [(0.6845 v U t logUM + (0.3423 log n- 5.277) v IT) |u | 


- 0.3423 nv U 4 t ] 


+ (0.4641 n 2 U 4 + 5.277 v 2 U 2 ) |U I + 1 .9525 nv U 4 /( v uV 
v T T' ' T* T T I 


dU 

d5 


U - | [-0.5 it Cl 


[-e C5 (3.39 DELP + 5.1) + (2.5 log C6 + 5. I) 2 V 
L w 


+ 2.5 log C6 + 5 . 1 ]/ 6 


+ 5.0 


V w (2.5 log C7 5.1) 
v C7 


°* 37 -C8 

+ L (3.39 ETAP + 1 .71 )e 

v 


2.5 U 


V C7 




0.5 nU 


— cr " c5 


e~^° U [1.2543 Cl U 3 6 4 
T l T 


- (1.2543 nuU 3 + 3.6387 Cl vU 2 )fi 3 
T T 


+ e C5 [Cl v 2 U (51.51 V + 7.6) 6 2 

L T W 



109 


+ [C2 C3 (-26.01 V + 5.1) + Cl y 3 (26.01 V +5.1) 

1 w w 

- 25.5 C3 V - 2.5 C3]6 

w 1 

+ C2 C4 (-26.01 V -5.1)1 
w 1 

+ e~ C5 log C6 [ci v IM38.0 V w + 2.5) fi 2 

+ C2 C3 (-25.5 V - 2.5) + Cl V 3 (25.5 V + 2.5) - 12.5 C3 V 1<5 
w w w J 

+ C2 C4 (-25.5 V - 2.5) ] 
w J 

+ e C ^ log 2 C6 [6,25 Cl v 2 V w fi 2 

+ (6.25 Cl v 3 V - 6.25 C2 C3 V )6 - 6.25 C2 C4 V 1 
w w w J 

+ [3.39 nv IT U 2 (C2 + 0.2487) - 9.993 Cl v 2 U.J6 2 

+ [8.49 C2 C3 + 1 .503 C3 - 5.1 Cl v 3 ]« 

2 2 3 

+ 5.1 C2 C4 ]/{ 2 v ttU t 6 + 2v ir6) 

- 0.5 O e (6 2 Cl - Cl nir6)/ir6 2 | 



no 


THIRD = 


FOURTH *= 


>.5 <U ^(-0. 5 it Cl 


[(2.5 log C6 + 5 . 1 ) V + 2.5 log C6 


(3.39 DELP + 5.1}e _C5 + 5.1 )/6 


5.0 (2.5 log C7 + 5.1 )U V 
T w 

C7 v 


0.37 U 

- 1 (3.39 ETAP + 5 . 1 )e _C 


3.39 U 


T -C8 
e 


• 5 U 0.5 it U ) 
+ £ 


v C7 


6 Cl - mr „ dU 
6 j it dU 




5.0 U V 6(2.5 log C6 + 5.1) 
U t [.0.5(C2 - »][■ T M V | U ; |- C6 1 


_C5 

+ 0.37 U^ 6 (3.39 DELP + 5.1)e /(v |u | ) 


- 3.39 C9 6 e" C5 + 2.5 

C6 


5.0 C9 n(2.5 log C7 + 5.1 )V 
v 

~ C7 


+ 0.37 C9 n(3.39 ETAP + 5.1 )e 


-C8 


- „ __ -C8 2.5 n C9 

- 3.39 C9 n e + ^ 



1 1 1 


2 

- 0.5(1 - C2)[v w (2.5 log C6 + 5.1) 

-C5 , 

+ 2.5 log C6 - (3.39 C6 + 5.1 )e + 5.1 ] 

2 — C8 

+ (2.5 log C7 + 5.1) V + 2.5 C7 - (3.39 ETAP + 5.1 )e +5.1 

w 


FIFTH = U f0.5 nn Cl fv (2.5 log C6 + 5.1 ) 2 
T 1 1 W 

-C8 2 

+ 2.5 log C6 - (3.39 DELP + 5.1 )e + 5.1 )/ 6 

- 0.5(1 - C2 ) ] [5.0 |uj ( 2 . 5 log C6 + 5.1 )/( v C6) 

-C8 2 * 5 l U J 0,5 nilU e 01 

+ 0.37 |U I (3.39 DELP + 4.1 )e /v+ — 1 

T V Cb J i 

SIXTH = 0.5 (1 - C2) 


where 


ETAP 


U n/ v 

T 


DELP * U 

T 


6/ v 


Cl 



^itn-v 

C2 = cos(-^J 



1 12 


„ 2 
C3 = nv tiu 


C4 = n TT 


C5 = 0.37 DELP 


C6 = DELP + 1 


C7 = ETAP + 1 


C8 = 0.37 ETAP 


U 




End of Document 



